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The mechanical response of static, unconfined, overcompressed face centred cubic, granular arrays 
is studied using large-scale, discrete element method simulations. Specifically, the stress response 
due to the application of a localised force perturbation - the Green function technique - is ob¬ 
tained in granular packings generated over several orders of magnitude in both the particle friction 
coefficient and the applied forcing. We observe crossover behaviour in the mechanical state of the 
system characterised by the changing nature of the resulting stress response. The transition between 
anisotropic and isotropic stress response exhibits critical-like features through the identification of 
a diverging length scale that distinguishes the spatial extent of anisotropic regions from those that 
display isotropic behaviour. A multidimensional phase diagram is constructed that parameterises 
the response of the system due to changing friction and force perturbations. 


I. INTRODUCTION 

Granular materials are ubiquitous throughout indus¬ 
try and nature and constitute myriad materials rang¬ 
ing from snowpacks, planetary regoliths, silos filled with 
grain, through to pharmaceutical pills and ceramic com¬ 
ponents. Yet predicting the mechanical properties of a 
given granular material, such as a packing of frictional 
grains, is difficult due to the many factors that influ¬ 
ence the stability of the packing. These can include 
microscopic features such as particle roughness, shape, 
and softness, to more macroscopic issues like the overall 
structural arrangement of the grains, density or packing 
fraction of the granular medium, and the properties of 
the container holding the material. Indeed, through a 
number of experimental and theoretical and com¬ 

puter simulation studies EHm, it has emerged that a 
granular medium can appear to behave not only as a me¬ 
chanically stable continuous elastic solid but can also ex¬ 
hibit strongly non-linear features characterised by highly 
heterogeneous stress properties. Though more recently, 
simulations that take into account the full preparation 
history of pouring a granular piling under gravity find 
that traditional elastoplastic continuum models are able 
to reproduce the stress-dip phenomena [12]. However, in 
the absence of such knowledge it is unclear how best to 
categorise any given granular packing. Yet, despite these 
complications, there are several methods that can be used 
to probe the mechanical robustness of the material. 

One protocol readily implemented that provides just 
such information on the nature of the stress state of a 
granular medium is the Green function technique: apply 
a localised, point force to the packing and observe the re¬ 
sulting stress profile in response to the perturbation m- 
A schematic of the procedure is shown in Figure This 
technique has been implemented in experiments [2] and 
computer simulations igiiiliii] to show that two dimen- 
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sional (2d) disc packings tend to agree with traditional, 
isotropic, elasticity theory nnniniin] when particle fric¬ 
tion is large. That is to say, the maximum in the mea¬ 
sured stress profile remains directly ‘below’ the location 
of the perturbing force and is often termed “one-peak” 
response in reference to two dimensional studies. This 
becomes a single lobe in three dimensions (3d). How¬ 
ever, for smoother particles, the stress response becomes 
highly anisotropic, exhibiting a “twin-peak” (ring profile 
in three dimensional packings) of maximum stress that is 
not directly in line with the applied force. Nevertheless, 
despite these differences, either type of response is gener¬ 
ally regarded to belong to solutions of the elastic-elliptic 
class [HiiiaiiH]- Furthermore, some studies [21 in] have 
also hinted at the possibility that the response within a 
grain pile may change character as one looks further from 
the location of the perturbing source. In another study 
m, a transition in the stress response within granular 
arrays is associated with a length scale though there has 
not been an explicit identification of the length scale nor 
the dependence on friction coefficient. This has impor¬ 
tant consequences as to whether a given material will in¬ 
deed conform to a more isotropic or anisotropic picture, 
and raises the possibility that the stress state of a partic¬ 
ulate medium may be tailored to respond differently on 
varying length scales. As yet there remains no method 
to predict a priori what the expected stress response of a 
three dimensional granular packing is and how the stress 
state of the system varies within the packing itself. 

It is precisely this issue that we address here. In an 
effort to reduce the large parameter space of all possible 
influences, we focus on what we reason are the two most 
dominant effects: particle friction and magnitude of the 
perturbing force. (The role of structural disorder was ad¬ 
dressed previously [20].) This study is needed to provide 
a guide in the design of granular solids where they can 
behave as a stable solid or exhibit strongly non-linear 
features with tailored mechanical properties. From this 
study, length scales distinguishing between anisotropic 
from isotropic behaviours can be estimated. The organ¬ 
isation of this paper is as follows. The Boussinesq equa- 
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FIG. 1. Implementing the force perturbation protocol referred 
to as the Green function technique. In this illustration, the 
arrow represents a downward, vertical point force applied at 
the middle on the top of a slab of material under interrogation. 
Possible stress response profiles are represented by the arcs 
and directed lines. The mechanical state of the systems is 
then inferred through monitoring the resulting stress profile 
along the a:-direction at different slices in z. Note, this is 
merely a 2d visualisation of the 3d problem we discuss here. 
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The crucial feature of the Green function approach is 
in reducing the problem to a material parameter free 
equation, Eq. The procedure we follow is through di¬ 
rect comparison of our simulation data with those of the 
Boussinesq result, Eq. When we have “single-peak” 
behaviour (in three dimensions, a lobed profile) we at¬ 
tempt to fit our results. We repeat this process over our 
range of parameters in friction coefficients and force per¬ 
turbations. 


III. SIMULATION PROCEDURE 


tions are briefly discussed in section |n| The manner in 
which we implement the packing preparation method and 


the Green function technique are described in section III 
and summarised in Ref. [21]. The results of this paper 
are presented in section [TVl followed by closing discussions 


and conclusions, section W 


II. BOUSSINESQ EQUATIONS 


Isotropic linear elastic theory is based on the solutions 
of the elliptic differential equations of mechanical equilib¬ 
rium m- Eor a point force (Eapp) applied in downward 
z direction to an elastic half space - the Green function 
technique (see Eigure - the solutions of the elliptic 
mechanical equilibrium differential equations provide the 
displacement field 
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where, Ux,y,z are the x, and 2 ; displacement field com¬ 
ponents and r = -h 2 :^. and G are bulk mate¬ 

rial parameters, the Poisson ratio and modulus of rigidity, 
respectively, u is restricted to the values between — 1 and 
0.5 HI]. The above equations, Eqs.[^ express a displace¬ 
ment vector at an arbitrarily chosen point inside a solid 
in response to the localised force perturbation Fapp. 

By implementing the displacement field components 
into the definitions of strain and applying Hooke’s law 
between stress and strain, we arrive at the stress response 
components known as Boussinesq equations (in three di¬ 
mensions) [ini [T7| . The origin of the coordinate system 
is at the application point of the applied forcing. Since 
the force is in the z direction the stress component char¬ 
acterising the response of a system is azz^ which we refer 


A. Packing Generation Protocol 


We employ simulation domains with periodic bound¬ 
ary conditions to reduce the effects of confining walls 
which can strongly influence the nature of the stress re¬ 
sponse in small systems [21]. Erom our previous studies 
on confined systems, we found that for a granular pack¬ 
ing that exhibits an isotropic, elastic-like response, the re¬ 
sponse behaviour converged to the theoretical predictions 
for systems of linear dimension, L > 20(i, for particle di¬ 
ameter d. Therefore, for this study we generated particle 
configurations of size L ^ 60d, as a system that is suitable 
to mimic the response within a semi-infinite medium. We 
placed N = 256000, monodisperse, inelastic spheres at 
the lattice points of a face centred cubic (ECC) array in 
the absence of gravity. The packing was then isotropi¬ 
cally compressed under uniform pressure until the system 
reached a packing fraction, (j) = 0.742, which is slightly 
above the hard-sphere ECC value, (j)FCC = 7r/\/l8- Par¬ 
ticles interact only on contact through a linear force law 
characterised by stiffness parameters, kn and kt^ for inter¬ 
actions normal and tangential to the contact plane where 
the maximum friction force is determined by the particle 
friction coefficient fi. Eor this study, the particle Poisson 
ratio was set to, kt/k^ = 1. Eriction is included using 
a standard granular dynamics routine [22]. The normal 
and tangential contact forces between particles (E^,t) are 
given by. 


Fji — kji^ij rUeffy "^nij 

Ft = -ktSij - 


(3) 


where Aij = {rij — d) is the surface compression of 
two contacting particles and Vn,t^j are the relative 
normal and tangential velocities of two particles. Vij 
is the relative distance between particles i at position 
Vi and j at position Vj. The effective mass, rueff is 
equal to m/2, for particle mass m. are the normal 
and tangential damping factors, which together with 
the kn,t and m set the coefficient of restitution of 
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the particles, which here is fixed to ^ 0 [23]. Sij 
represents the integrated displacement of the contact 
point of the particles i and j while the two particles 
remain in contact. The Coulomb condition is written 
as Ft < fiFn- We generated static packings for different 
friction, ranging over several orders of magnitude, /r = 
0,0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.5,1,10. 
Below, we present our results in simulation units where 
mass, length, and force values are appropriately scaled 
by the particle properties: mass m = 1, diameter d = 1, 
and stiffness k = = kt = 1. 


B. Force Perturbation Protocol 
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FIG. 2. Simulation snapshots of the force perturbation in 
systems that exhibit, (a) an anisotropic, ring-like (twin-peak) 
response /i = 0 and (b) an isotropic, single lobe (single-peak) 
response, ji — 1. Each panel focuses on the central region 
of the packings with particles shaded by their displacement 
relative to the unperturbed initial state. A force perturba¬ 
tion was applied at the centre of the box as described in the 
text. Shading represents magnitude of the displacement (rel¬ 
ative to the average for that system): blue end of spectrum 
corresponds to larger displacements. For the FCC particle ar¬ 
rangement at an overcompressed packing fraction 0 = 0.742 
used in this study. 


The localised force perturbations were implemented 
in the following manner. We first located an approxi¬ 
mately spherical cluster of particles within the central 
region of the packing that typically consisted of < 0.3% 
of the total number of particles in the system. We 
then increased their diameters by a specified amount 
and allowed the packing to relax back into a mechan¬ 
ically stable state keeping the same periodic boundary 
conditions. The effect of this local swelling of particles 
is equivalent to applying a localised and isotropic per¬ 
turbing force, Tapp, at the centre of the packing, that 
we varied over several orders in magnitude: Tapp = 
0.0005,0.001,0.0025,0.005,0.01,0.025,0.05,0.1,0.25 kd. 
This procedure, initially implemented in Ref.[20l was mo¬ 
tivated by several considerations. Firstly, given that we 
exclude gravity in this work, we preferred to avoid im¬ 
posing a uni-directional force that would break the sym¬ 
metry of our set up. We also reasoned that this protocol 
offers a practical means to implement force perturbations 
in real systems. As a case in point, during the course of 
this work, a recent experimental work has implemented 


a similar procedure in studies of the mechanical proper¬ 
ties of two-dimensional granular packings [24]. We have 
verified that the method of increasing the size of parti¬ 
cles, located within a small region centred at origin, does 
indeed result in a corresponding force magnitude equal 
to Tapp, by measuring the average force over a spherical 
shell in the vicinity of the inflated region. Representative 
snapshots of how the force perturbation influences the 
packing are shown in Figure [^ These images display the 
particles occupying a sample, cubic central portion of the 
packing (out to lOd) for /i = 0 and Tapp = 0.25 (Figure 
l^a)) and /i = 1 and Tapp = 0.02 (Figure [^b)). These 
examples were chosen to clearly illustrate the change 
in response. The magnitude of particle displacements 
relative to their initial, unperturbed positions are indi¬ 
cated by a shading scale, with the blue end of the spec¬ 
trum indicating larger displacements. These microdis¬ 
placements of the particles in response to the applied 
force indicate the reasons behind the differences between 
the anisotropic and isotropic response functions. In the 
anisotropic case, maximal particle displacements are di¬ 
rected along an open cone, whereas in the isotropic case 
the displacements ultimately form a closed ‘ball’ ema¬ 
nating away from the region of the perturbing source. 
We have verified that for the system size studied here, 
our results match exactly the predictions of Eq. when 
we observe an isotropic response (which occurs for large 
friction and small applied force - see Figure |^. 

We define the stress response as the difference between 
the final and initial stress states. We computed the stress 
in the packing before, di, and after, df, the imposition of 
the perturbation, so that the resulting stress response is, 

cr = (7f - (Tj. (4) 

C. Coarse Graining Procedure 

To be able to accurately compare our simulation stress 
response with elasticity theories, we implement a stress 
computation method that connects our microscopic, dis¬ 
cretized calculations with macroscopic, continuum length 
scales. A stress expression with a coarse-graining pro¬ 
cedure developed by Goldhirsch and colleagues [mus] 
gives continuum stress components. We chose a Gaus¬ 
sian coarse-graining function, so that the components 
(=x,^, 2 :), of the stress tensor cr(r), are obtained 

from, 

(5) 

where uj is the coarse-graining length scale and are 
the contact forces between particles i and j. This ex¬ 
pression is normalised. The value of the coarse-graining 
length scale, cj, is important in terms of obtaining a fine 
resolution and yet continuum description of stress compo¬ 
nents. Here, the particle diameter, d, is a suitable choice. 
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The coarse grained stress expression, Eq. is calculated 
on a spatial grid discretized by grid size of 0.2d. 


IV. RESULTS 

A. Stress Maps 

To make our results connect visually with existing ex¬ 
periments, we present our results as stress maps in dif¬ 
ferent planes with respect to the perturbation source. 
We also plot the circularly-averaged stress as two dimen¬ 
sional profiles across the central section from planes that 
are perpendicular to the direction of the force perturba¬ 
tion. In Figurej^we show stress profiles for packings with 
different friction subject to an applied force perturbation 
of a relatively large, fixed magnitude Fapp = 0.25 kd. 
The panels of Figure are arranged as follows: each row 
corresponds to a single friction coefficient while each col¬ 
umn shows data at different distances from the perturb¬ 
ing source. We refer to these planes as Top’, ‘middle’, 
and ‘bottom’ referring to how far from the perturbation 
source we view the stress, with the top layer closest to 
the source and bottom furthest away. The left column 
is stress profile closest to the source (‘top’) and the right 
column furthest away (‘bottom’). Immediately we see 
how friction plays a major role in determining the re¬ 
sponse of the system. Consistent with previous studies 
[9] we find that friction enhances the elastic nature of 
the response whereby the response profile remains “sin¬ 
gle peaked” at all distances from the source, with linear, 
isotropic elastic broadening of the profile. At zero fric¬ 
tion, the profile is strongly anisotropic away from the 
source with a ring-like or “twin-peak” response, reflect¬ 
ing the lattice structure of the underlying packing, that 
persists at the largest distances. However, for interme¬ 
diate friction values, the response initially resembles the 
zero-friction state out to a finite distance from the source 
then transitions to a more isotropic, single-peak response 
at larger distances. 

B. Stress Profiles 

To provide a more quantitative comparison between 
the different plots shown in Figure as well as to present 
our results in a manner that is closer in spirit to similar 
data obtained from 2d experiments, we show one dimen¬ 
sional stress profiles in Figures and[^ These stress 
profiles are obtained by circularly averaging the ‘bottom’ 
stress maps. Figure shows the stress response profile, 
cr, at the bottom of the box for a relatively large, applied 
force Fapp = 0.2bkd, for packings with different friction 
coefficients. These profiles clearly indicate the role of fric¬ 
tion in determining the nature of the stress response. Low 
friction results in a twin-peak profile of the anisotropic 
elastic variety, that transitions to a single-peak profile 
at larger friction corresponding to the isotropic, linear 


elastic Boussinesq theory m- These results are consis¬ 
tent with similar studies performed in two dimensional 
geometries [9]. 

For a given friction coefficient the stress response is 
also sensitive to the magnitude of the perturbing force. 
We illustrate this point in Figure where we show the 
stress profiles at the bottom of a packing with /i = 0.01, 
for different values of Fapp. Again, the common theme 
here is that the response changes character depending on 
Fapp. At smaller forces, there is little or no rearrange¬ 
ment in the particle contacts, and the underlying force 
network remains intact leading to a single peak profile. 
Whereas, for larger forces, a small fraction of contacts 
are broken, thereby causing the stresses to be propa¬ 
gated along better defined paths, leading to twin peak 
profiles. For the case where a single-peak response is ob¬ 
tained, we check that these profiles are consistent with 
isotropic, linear elastic theory, by fitting the profiles to 
the Boussinesq equation, Eq. as shown in Figure 
and extract the full width at half height. We note that 
there are no free parameters in this fitting procedure as 
the applied forcing is given as an input. As shown in 
Figure we display our analysis as the ratio of the fit¬ 
ting to the expected theoretical result. Thus we confirm 
that the systems studied here do behave as an isotropic, 
elastic continuous medium over a range of parameters in 
friction and applied forcing. 

To gain better insight into this transitional stress re¬ 
sponse we show the stress profiles in the plane parallel 
to the direction of the initial perturbation in Figure 
These two-dimensional plots were obtained by spheri¬ 
cally averaging the stress response about the location of 
the perturbation point. This view more clearly demon¬ 
strates how an initial anisotropic response - characterised 
by the twin-peak stress profile - transitions towards one 
with more traditional isotropic features - characterised 
by single lobe response profile - and ultimately how the 
extent of the force perturbation is eventually, in some 
sense, dissipated into the background stress. Closer in¬ 
spection of the stress response, as visualised in Figure 
reveals that anisotropic response features can extend ar¬ 
bitrarily far into the packing depending on the packing 
parameters. 


C. Response Scaling 

To characterise the transition between anisotropic and 
isotropic elastic stress response behaviours we designed 
the following procedure that identified how the maximum 
stress transmitted through the system deviated from tra¬ 
ditional isotropic behaviour. For an isotropic, linear 
elastic continuum described by the Boussinesq equation, 
the maximum stress occurs in a line directly “below”, 
the point of force perturbation. Our circularly averaged 
stress analysis allows us to scan over angle 6>, to identify 
the direction of maximal stress transmission and mon¬ 
itor this response out to a radial distance, r, from the 
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FIG. 3. The stress response viewed in planes normal to the direction of the perturbing force. The origin of the coordinate 
system is defined as the location of the perturbing force. Here, the horizontal slices are denoted xy-planes and the different 
columns correspond to different slices in the z direction. The three rows correspond to different particle friction: /x = 0 (top), 
/i = 0.05 (middle), and /x = 10 (bottom). Within each row, columns show the stress at different distances from the location of 
the perturbing force: left column is closest to the force. The shading scale is the same for each row but different among the 
columns to better reveal the stress profiles. In all cases, the applied force magnitude, Fapp = ^.2bkd. 



FIG. 4. Gircularly averaged stress profiles a, at the ‘bottom’ 
layer due to applied force iTpp = 0.25/cd, for different friction 
coefficient /x, as indicated in key. 


FIG. 5. Gircularly averaged stress profiles a, at the bottom 
layer due to forces, at /x = 0.01, for different applied force 

p 

-L app' 


perturbation source. This construct is shown in Figure 

The anisotropic stress response in our ‘^D geometry 
transmits through a cone-like structure. Utilising the az¬ 
imuthal symmetry of this geometry we angularly average 
the stress response, over angular elements A6> = 4°, lead¬ 
ing to a 2D projection of the result as shown in Figure 


a precise measure of the direction of maximum stress re¬ 


10 The two vertical lines in Figure pT] provides us with 


spouse relative to the direction for the expected isotropic 
maximum stress direction (6> = 0). We define the angular 
direction 0^ of maximal stress response as the average of 
these two peak locations. To identify a length scale that 
quantifies the nature of the mechanical state we measure 
the distance, r, from the origin of the localised force, in 
the direction of maximal stress, until that value of the 
stress is subsumed into the background stress. Starting 
from the angular slice closest to the perturbation source 
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FIG. 6. The influence of friction on the stress response. Profiles are viewed in planes parallel to the direction of the applied 
force perturbation, for a force magnitude Fapp = 0.25kd. Panels are particle friction values: (top left) /x = 0, (top middle) 
/i = 0.01, (top right) /i = 0.05, (bottom left) fi = 0.5, (bottom middle) /x = 1, and (bottom right) /x = 10. 




FIG. 7. Gircularly averaged stress profiles cr, at the bottom 
layer due to forces, with ji — 10, for different applied force 
Fapp (color dashed lines) and the fittings to the Boussinesq 
equation (solid black lines). 


we scan elemental layers, layer-by-layer, away from the 
source to track the direction of the maximum stress. At 
each layer we compare the value of the stress in this di¬ 
rection with the stress in the same radial layer in the 
direction below the perturbing source, 6> = 0. Thus, we 
define the characteristic length scale as the distance at 
which these two stress values become indistinguishable, 
as indicated by the vertical line in Figure EH 

For the extensive parameter space in /x and Fapp ex¬ 
plored here, we find that both the characteristic length 
scale and angle exhibit power law dependence on 
friction coefficient of the form 

S,H oc (6) 

9^ cx for /i > 0 (7) 

with the power-law exponents, a = 0.5 ± 0.06 and 
(3 = 0.37T0.07, over the range of friction and force magni¬ 
tudes studied. Thus, the anisotropic-isotropic transition 


FIG. 8. Full width at half height. A, of the stress profiles, for 
a range of applied force perturbations and friction coefficients, 
when a single-peak stress response is observed. Ab — 1.113z 
is the theoretical result. For each friction coefficient, data is 
shifted by one unit for convenience, (/x = 10 at top descending 
in order to /x = 0.02 at bottom.) For /x = 0.05 and Fapp — 
0.25A;d and fi — 0.02 and Fapp = O.lkd, larger full widths 
of the profiles indicating a deviation from the theory. For 
/X = 0.02 and Fapp = 0.25A;d, the profile is shown in Figure]^ 
and appears as a twin-peak profile. 


for stress response of a granular array is characterised 
by critical-like features accompanied by a length scale 
that grows as friction decreases. A subset of the data 
(for visual clarity) is presented in Figure 12 The very 
presence of a diverging length scale that distinguishes 
between different mechanical phases suggests that the 
anisotropic-isotropic transition for stress response of a 
granular array is characterised by critical-like features 
that are controlled by the friction coefficient. Curiously, 
this behaviour mimics the manner in which a thermo¬ 
dynamic system behaves in the vicinity of a continuous 
phase transition where now the friction coefficient plays 
the role of a reduced temperature. In this sense, the 
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FIG. 9. Scanning the circularly average stress response over 
angle 0 at radial slices r, to identify the direction of the 
anisotropic response profile relative to the direction for the 
isotropic case (dashed line). 



FIG. 10. Gircularly averaged stress profile, a, for a system 
exhibiting an anisotropic, twin-peak prohle (Fkpp = 0.25/cd 
and /i = 0.05). The profile peaks, identified with vertical 
lines, indicate the directions of the maximal stress and are 
denoted 0^ in the text. 



FIG. 11. The circularly averaged stress response profile cr, 
in the direction 0^ of maximal stress (red solid line) and the 
background stress (green dashed line) in the direction ^ = 0, 
identified in Figure|^(for, iTpp = 0.2bkd and — 0.05). This 
representative dataset shows how the characteristic length 
scale, at which the directional stress crosses the back¬ 
ground stress is identified by the vertical line. 


length scale defined here measures the spatial range of 
fluctuations in the stress over and above the background 
isotropic stress. Thus, to some extent the anisotropic- 
isotropic transition resembles a one-dimensional Ising 
model, where the zero-friction limit corresponds to the 
ordered phase and infinite friction is the analogue of the 
disordered phase, stable fixed point. For any friction 
‘perturbations’, the response is readily driven away from 
the zero-friction, anisotropic phase, and always ends up 
exhibiting isotropic stress response for large enough sys¬ 
tem size. 

Upon collating the data over our entire parameter 
space: fi = {0,10.0} and Fapp = {0.0005,0.25}, we 
are able to map the stress response behaviour onto the 
pseudo-phase diagrams shown in Figure These di¬ 
agrams delineate regions of parameter space where an 
anisotropic response is expected (region below the bound¬ 
ary) and those parameters which are more likely to result 
in an isotropic, linear elastic response (above boundary). 
Anisotropic behaviour emerges for low friction and high 
forcing where the stress peaks extend far into the pack¬ 
ing and are directed away from perturbation direction. 
It is only in the /i = 0 limit that the anisotropic length 
scale spans the simulation domain. Conversely, isotropic 
character is recovered for mildly frictional materials and 
moderate forcing. 

Combining these two plots we present the anisotropic- 
isotropic transition in stress response as parameterised 
by the effective length scale, = ^^sin6>, where is 
normalised by the simulation box size, in Figure [T4| The 
boundary distinguishes between the stress state within a 
single system. For points below the boundary, the gran¬ 
ular medium behaves as an anisotropic, elastic material, 
where stresses are transmitted through the system along 
directed pathways commensurate with the structure of 
the packing. Above the boundary, the material will re¬ 
sort to behaviour more consistent with an isotropic, elas¬ 
tic solid. Thus, small systems subject to relatively large 
forces will appear to behave as anisotropic materials that 
may result in failure along the directions of the largest 
stress propagation. Whereas large frictional aggregates 
can be considered, as is well-known from everyday expe¬ 
rience tells us, as isotropic elastic media. 


V. CONCLUSIONS 

Large-scale computer simulations of overcompressed 
FCC granular arrays over a wider range in friction co¬ 
efficient and point force perturbations were carried out 
to study the stress response inside these arrays. Our re¬ 
sults indicate that packings with large particle friction 
subject to small localised forcings, exhibit a single lobe 
stress response consistent with classical, isotropic elastic¬ 
ity. Whereas, packings with low friction and subject to 
large localised force magnitudes exhibit ring-like stress 
responses, a hallmark of anisotropic elastic behaviour. 
However, this anisotropic state diminishes with distance 















FIG. 12. Influence of friction and force magnitude on the transition between isotropic and anisotropic stress response. (Left) 
Length scale characterising the distance from the location of the perturbing force at which the response crosses over from 
anisotropic to isotropic. The lines are power-law fits: ~ with a — 0.5 ± 0.06. The scaled length data points are 

collapsed onto the line in inset. (Right) Angle direction of the anisotropic stress transmission relative to the central, isotropic 
direction. The lines are power-law fits: 0^ ~ for /x > 0, with (3 — 0.37 ± 0.07. The scaled angles apart from uppermost 
data points are collapsed onto the line in inset. In both panels, the quoted force perturbations are in units of kd. 



FIG. 13. Phase diagrams of force perturbations, in applied force iTpp, and friction coefficient /x, for the (a) and (b) stress 
transmission direction 0. These 3D plots delineate between regions of isotropic, single-peak response (above) and anisotropic, 
twin-peak, response (below the dividing surface). 



FIG. 14. Three dimensional phase diagram, in the space of 
applied force Fapp and particle friction coefficient /x, for the 
anisotropic-isotropic transition of mechanical response within 
a granular array, parameterised by the effective length scale, 
sin^ where is normalised by simulation box size. 


from the perturbation source, within the same packing, 
and transitions towards a more elastic-like character at a 
length scale that depends on the friction coefficient. This 
characteristic length scale grows with decreasing friction. 
Ultimately, the limit of purely frictionless packings ap¬ 
pear to present as a pathological case in this problem 
- showing anisotropic behaviour out to the largest dis¬ 
tances. Any non-zero friction will ultimately send the 
system towards an isotropic elastic regime in the ther¬ 
modynamic limit. 

While the FCC arrays used in this study clearly en¬ 
hance the possibility of observing strongly anisotropic be¬ 
haviour, the presence of disorder is likely to suppress such 
features [9l[20]. Thus, we expect that a combination of 
varying structural disorder and friction provide the nec¬ 
essary ingredients for packings that can be constructed 
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with tailored mechanical properties. It might also be 
that as a disordered granular system is brought closer to 
its point of marginal rigidity, it may become increasingly 
fragile [261127], and remain as an anisotropic medium at 
all friction values, thus providing a broad range of possi¬ 


ble mechanical states to choose from. 
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